#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include "gmp.h"

//Parameters for GMP lib
#define ITER            25
#define PRECISION       100000000
#define BASE            10
#define DIGITS          30000001
//Number of threads
#define N_THREADS       4
//Number of mutexes
#define N_MULTEXES      3

//Global variables of each iteration of variables a,b,t,p
mpf_t a0[ITER + 1], b0[ITER + 1], t0[ITER + 1], p0[ITER + 1];
//Iteration set by user
int iter;

//Declaration of mutexes
pthread_mutex_t mutex[N_MULTEXES];

//Vector used to control the actual iteration
//This vector is used to know in what iteration each variable is
int iteration[4] = {0};

//Calculates an1 = ((a0 + b0) / 2.0);
void* an1();

//Calculates bn1 = sqrt(a0 * b0);
void* bn1();

//Calculates tn1 = tn - pn(pow(an-an1,2))
void* tn1();

//Calculates  pn1 =  2.0 * p0
void* pn1();

//Main function
int main(int argc, char *argv[]) {

  //Verify number of iteration
    if(argc == 2)
        iter = iter < 25 ? atoi(argv[1]) : ITER;
    else iter = ITER;

    fprintf(stderr, "Calculando Gauss-Legendre Paralelo com %d iteracoes ...\n",iter);

    //Create threads
    pthread_t threads[N_THREADS];
    //Aux variables to help in the final calculus
    mpf_t tmp, aux, result;
    // Set default precison
    mpf_set_default_prec(PRECISION);
    //Initialize a, b, t and p vectors
    int i = 0;
    for (i = 0; i <= iter; i++)
        mpf_inits(a0[i], b0[i], t0[i], p0[i], NULL);
    //Initialize aux variables
    mpf_inits(tmp, aux, result, NULL);

    //Set a0 initial valor
    //    a0 = 1.0;
    mpf_init_set_ui(a0[0], 1);

    //Set b0 initial valor
    //    b0 = 1.0 / M_SQRT2;
    mpf_init_set_ui(b0[0], 1);
    mpf_init_set_ui(b0[1], 2);
    mpf_sqrt(b0[1], b0[1]);
    mpf_div(b0[0], b0[0], b0[1]);
    //Set t0 initial valor
    //    t0 = 1.0 / 4.0;
    mpf_init_set_ui(t0[0], 1);
    mpf_init_set_ui(t0[1], 4);
    mpf_div(t0[0], t0[0], t0[1]);


    //Set p0 initial valor
    //    p0 = 1.0;
    mpf_init_set_ui(p0[0], 1);

    //Initialize mutexes
    for (i = 0; i < N_MULTEXES; i++)
        pthread_mutex_init(&mutex[i], NULL);

    //Initialize every thread
    pthread_create(&threads[0], NULL, an1, NULL);
    pthread_create(&threads[1], NULL, bn1, NULL);
    pthread_create(&threads[2], NULL, tn1, NULL);
    pthread_create(&threads[3], NULL, pn1, NULL);

    //Wait for the end of every thread
    for (i = 0; i < N_THREADS; i++)
        pthread_join(threads[i], NULL);


    //Destroy every mutexes
    for (i = 0; i < N_MULTEXES; i++)
        pthread_mutex_destroy(&mutex[i]);


    //Calculus of pi = pow(an+bn,2)/(4*tn)
    mpf_add(tmp, a0[iter - 1], b0[iter - 1]);
    mpf_pow_ui(tmp, tmp, 2);


    mpf_mul_ui(aux, t0[iter - 1], 4);
    mpf_div(result, tmp, aux);

    //Print the result
    mpf_out_str(stdout, BASE, DIGITS, result);

    //Clear the vectors
    for (i = 0; i <= iter; i++)
        mpf_clears(a0[i], b0[i], t0[i], p0[i], NULL);
    //Clear the aux variables
    mpf_clears(result, tmp, aux, NULL);


    return 0;
}


// an1 =  ((a0 + b0) / 2.0);
void* an1() {
    //Temporary variables to help
    int i, baux = 0;

    for (i = 0; i < iter; i++) {
        //This while is to verify in what iteration b is
        //While b iteration is lower than a iteration, a waits for b
        while (baux < i) {
            pthread_mutex_lock(&mutex[1]);
            baux = iteration[1];
            pthread_mutex_unlock(&mutex[1]);
        }
        //This is the calculus of a
        mpf_add(a0[i + 1], a0[i], b0[i]);
        mpf_div_ui(a0[i + 1], a0[i + 1], 2);

        //This is for update the actual iteration of variable a
        pthread_mutex_lock(&mutex[0]);
        iteration[0] = i + 1;
        pthread_mutex_unlock(&mutex[0]);
    }
    return;
}

//bn1 = sqrt(a0 * b0);
void* bn1() {
    //Temporary variables to help
    int i = 0, aaux = 0;

    for (i = 0; i < iter; i++) {
        //This while is to verify in what iteration a is
        //While a iteration is lower than a iteration, b waits for a
        while (aaux < i) {
            pthread_mutex_lock(&mutex[0]);
            aaux = iteration[0];
            pthread_mutex_unlock(&mutex[0]);
        }

        //This is the calculus of a
        mpf_mul(b0[i + 1], a0[i], b0[i]);
        mpf_sqrt(b0[i + 1], b0[i + 1]);

        //This is for update the actual iteration of variable b
        pthread_mutex_lock(&mutex[1]);
        iteration[1] = i + 1;
        pthread_mutex_unlock(&mutex[1]);
    }
    return;
}

//tn1 = tn - pn(pow(an-an1,2))
void* tn1() {
    //Temporary variables to help
    int i, aaux = 0, paux = 0;

    for (i = 0; i < iter; i++) {
        //This while is to verify in what iteration a is
        //While an+1 iteration is lower than t iteration, t waits for an+1
        while (aaux <= i) {
            pthread_mutex_lock(&mutex[0]);
            aaux = iteration[0];
            pthread_mutex_unlock(&mutex[0]);
        }

        //This while is to verify in what iteration a is
        //While p iteration is lower than t iteration, t waits for p
        while (paux < i) {
            pthread_mutex_lock(&mutex[2]);
            paux = iteration[3];
            pthread_mutex_unlock(&mutex[2]);
        }


        ////This is the calculus of t
        //   tmp = pow(a0 - a1, 2);
        mpf_sub(t0[i + 1], a0[i], a0[i + 1]);
        mpf_pow_ui(t0[i + 1], t0[i + 1], 2);

        //    (t0 - (p0 * (tmp)))
        mpf_mul(t0[i + 1], p0[i], t0[i + 1]);
        mpf_sub(t0[i + 1], t0[i], t0[i + 1]);


        //This is for update the actual iteration of variable t
        iteration[2] = i + 1;

    }
    return;
}

//  pn1 =  2.0 * p0
void* pn1() {

    int i;
    for (i = 0; i < iter; i++) {
        mpf_mul_ui(p0[i + 1], p0[i], 2);

        //This is for update the actual iteration of variable p
        pthread_mutex_lock(&mutex[2]);
        iteration[3] = i + 1;
        pthread_mutex_unlock(&mutex[2]);
    }
    return;
}

